!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

! qdr_xxx - public functions
!
! 1D - function integration
module qdr_mod; 
use func_mod;
use vec_mod;
implicit none;

logical, save :: QUADR_TRACE = .TRUE.;

private;

! NOTE 9-Apr-2013: name any new public functions as qdr_xxx()
public  test_qdr_mod, qdr_step_p5 &
   & , qdr_max_orth_err, qdr_max_orth_err3, qdr_lcr_max_orth_err3 &
   & , qdr_intgrl, qdr_intgrl2, qdr_intgrl3 &
   & , qdr_lcr, qdr_intgrl_lcr, qdr_intgrl_lcr2, qdr_intgrl_lcr3&
   & , qdr_intgrl_lcr2DivR, qdr_intgrl_lcr2DivR2 ;
contains 
!___________________________________________________________
SUBROUTINE test_qdr_mod(); 
  character(LEN=64) :: name = "test_qdr_mod";
  call dbg_in(name);
  call test_qdr();
  call dbg_out(name); 
end SUBROUTINE
!___________________________________________________________
SUBROUTINE test_qdr();  
  integer :: nv;
  real(DP) :: minv, maxv, res, res2, MAX_DRV_ERR = 1.d-18;
  real(DP), pointer :: px(:);      
  real(DP), dimension(301) :: py, xq;      
  character(LEN=64) :: name = "test_qdr";
  call dbg_in(name);

  nv = 5;  minv = D0; maxv = D1;
  px => mem_vec(nv);  
  call grid_step(px, minv, maxv);
  call dbg_assert_zero("  5.px(1)=minv", px(1) - minv);
  call dbg_assert_zero("  5.px(2)=minv", px(2) - (maxv+minv)/4.);
  call dbg_assert_zero("  5.px(nv)=maxv", px(nv) - maxv);
  
  call qdr_step_p5(xq, px); call vec_dbg("xq=", xq); ! x-quadratures
  call dbg_assert_zero("  xq(1)", xq(1) - D05*7D0  /45D0);
  call dbg_assert_zero("  xq(2)", xq(2) - D05*32D0 /45D0);
  call dbg_assert_zero("  xq(3)", xq(3) - D05*12D0 /45D0);
  call dbg_assert_zero("  xq(4)", xq(4) - D05*32D0 /45D0);
  call dbg_assert_zero("  xq(5)", xq(5) - D05*7D0  /45D0);

  call fun_vec_polynom(py, px, (/D1, D0, D0/)); 
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("int 1.=", res, D1, EPS25);

  ! new x-grid-9
  call mem_del_vec(px);
  px => mem_vec(9);  
  call grid_step(px, D0, D1);
  call fun_vec_polynom(py, px, (/D1, D0, D0/)); 
  call qdr_step_p5(xq, px); ! x-quadratures
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("9.int 1.=", res, D1, EPS15);
  
  ! new x-grid-13
  call mem_del_vec(px);
  px => mem_vec(13);  
  call grid_step(px, D0, D1);
  call fun_vec_polynom(py, px, (/D1, D0, D0/)); 
  call qdr_step_p5(xq, px); ! x-quadratures
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("13.int 1.=", res, D1, EPS15);
  
  call fun_vec_polynom(py, px, (/D0, D1/)); 
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("13.x^2=", res, D05, EPS15);

  call fun_vec_polynom(py, px, (/D0, D0, D1/)); 
  res = qdr_intgrl(xq, py); 
  res2 = D1/D3;
  call dbg_assert_abs_diff_LE("13.x^3=", res, res2, EPS15);

  call fun_vec_polynom(py, px, (/D1, -D1, D1, -D1, D1/)); 
  res = qdr_intgrl(xq, py); 
  res2 = D1 - D05 + D1/D3 - D1/D4 + D1/D5;
  call dbg_assert_abs_diff_LE("13.1-1+1-1+1=", res, res2, EPS15);

  ! new x-grid-21
  call mem_del_vec(px);
  px => mem_vec(21);  
  call grid_step(px, D0, PI);
  call qdr_step_p5(xq, px); call vec_dbg("qdr_step_p5(xq, px)=", xq);
  
  !21-Mar-2013: =SIN(px) is not working.
  !py = SIN(px);  call vec_dbg("py=SIN(px)", py);
  call vec_sin(py, px);  call vec_dbg("py=SIN(px)", py);
  
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("21.sinx=", res, D2, EPS7);

  !py = COS(px); !NOTE: =COS(px) is not working.
  call vec_cos(py, px);  call vec_dbg("py=COS(px)", py);
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("21.(1-cosx)=", res, D0, EPS15);

  ! new x-grid-301
  call mem_del_vec(px);
  px => mem_vec(301);  
  call grid_step(px, D0, PI);
  call qdr_step_p5(xq, px); ! x-quadratures
  !py = SIN(px);
  call vec_sin(py, px);  call vec_dbg("py=SIN(px)", py);
  res = qdr_intgrl(xq, py); 
  call dbg_assert_abs_diff_LE("301.sinx=", res, D2, EPS14);

  py = D1-COS(px);  call vec_dbg("py=D1-COS(px)=", py);
  res = qdr_intgrl(xq, py);   call dbg_dp("res=", res);
  call dbg_assert_abs_diff_LE("21.(1-cosx)=", res, PI, EPS24);

  ! new x-grid-101
  call mem_del_vec(px);
  px => mem_vec(101);  
  call grid_step(px, -D1, D1);
  call qdr_step_p5(xq, px); ! x-quadratures
  !py = EXP(px);
  call vec_exp(py, px);  call vec_dbg("py=EXP(px)", py);
  res = qdr_intgrl(xq, py); 
  res2 = EXP(D1)-EXP(-D1);
  call dbg_assert_abs_diff_LE("101.exp=", res, res2, EPS12);

  call mem_del_vec(px);
end SUBROUTINE;
!___________________________________________________________
real(DP) function qdr_max_orth_err(xqdr, av) result(res);        
  real(DP), intent(IN) :: xqdr(:), av(:,:);     ! av - array of vectors  
  integer :: i, i2, min2, max2;
  real(DP) :: w;
  call dbg_str("-->qdr_max_orth_err()", QUADR_TRACE); 
  call mtrx_i2_range(min2, max2, av);
  res = D0;
  do i = min2, max2  
    call dbg_int("i=", i, QUADR_TRACE);
    do i2 = min2, i  
      w = qdr_intgrl2(xqdr, av(:, i), av(:, i2));
      if (i == i2)  w = w - D1;
      w = ABS(w);
      if (res < w)  then 
        res = w;  
        call dbg_int("i2=", i2, QUADR_TRACE);
        call dbg_dp("max_err=", res, QUADR_TRACE);
      endif 
    end do
  end do 
end function
!___________________________________________________________
!real(DP) function calc_max_orth_err3(xqdr, av3) result(res);        
!  real(DP), intent(IN) :: xqdr(:), av3(:,:,:);     ! av - array of vectors  
!  integer :: i3, min3, max3;
!  real(DP) :: w;
!  call dbg_str("-->qdr_max_orth_err()"); 
!  call mtrx_i2_range3(min3, max3, av3);
!  !call dbg_int("min3=", min3);
!  !call dbg_int("max3=", max3);
!  res = D0;
!  do i3 = min3, max3  
!    w = qdr_max_orth_err(xqdr, av3(:,:,i3)); 
!    !call dbg_dp("w=", w);
!    w = ABS(w);
!    if (res < w) then  
!      res = w;
!      call dbg_dp("new max err=", res);
!    endif 
!  end do 
!end function
!___________________________________________________________
real(DP) function qdr_max_orth_err3(xqdr, av3) result(res);        
  real(DP), intent(IN) :: xqdr(:);     ! av - array of vectors  
  type(ARR_VEC_PTR_T), pointer :: av3(:);  
  integer :: i3, min3, max3;
  real(DP) :: w;
  call dbg_str("-->calc_max_orth_err4()"); 
  call arr_vec_ptrs_range(min3, max3, av3);
  !call dbg_int("min3=", min3);
  !call dbg_int("max3=", max3);
  res = D0;
  do i3 = min3, max3  
    w = qdr_max_orth_err(xqdr, av3(i3)%av); 
    !call dbg_dp("w=", w);
    w = ABS(w);
    if (res < w) then  
      res = w;
      call dbg_dp("new max err=", res);
    endif 
  end do 
end function
!___________________________________________________________
real(DP) function qdr_lcr_max_orth_err3(lcr, av3) result(res);        
  type(GRID_LCR_T), intent(IN) :: lcr;
  type(ARR_VEC_PTR_T), pointer :: av3(:);  
  res = qdr_max_orth_err3(lcr%xqdr_cr2, av3);
end function
!___________________________________________________________
SUBROUTINE qdr_lcr(lcr);  
  type(GRID_LCR_T), intent(INOUT) :: lcr;
  call dbg_str("-->qdr_lcr()", QUADR_TRACE);
  call dbg_int("lcr%nx=", lcr%nx, QUADR_TRACE);
  lcr%xqdr => mem_vec(lcr%nx);
  call qdr_step_p5(lcr%xqdr, lcr%xv);
  lcr%xqdr_cr => mem_vec(lcr%nx);
  lcr%xqdr_cr = lcr%xqdr * lcr%cr;
  lcr%xqdr_cr2 => mem_vec(lcr%nx);
  lcr%xqdr_cr2 = lcr%xqdr_cr * lcr%cr;  
  lcr%xqdr_cr2DivR2 => mem_vec(lcr%nx);
  lcr%xqdr_cr2DivR2 = lcr%xqdr * lcr%crDivR * lcr%crDivR;  
  lcr%xqdr_cr2DivR => mem_vec(lcr%nx);
  lcr%xqdr_cr2DivR = lcr%xqdr_cr * lcr%crDivR;  
end SUBROUTINE
!___________________________________________________________
real(DP) function qdr_intgrl_lcr(lcr, fx) result(res);   
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP), intent(IN) :: fx(:);      
  res = qdr_intgrl(lcr%xqdr_cr2, fx);
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl_lcr2(lcr, fx, fx2) result(res);   
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP), intent(IN) :: fx(:), fx2(:);      
  res = qdr_intgrl2(lcr%xqdr_cr2, fx, fx2);
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl_lcr2DivR2(lcr, fx, fx2) result(res);   
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP), intent(IN) :: fx(:), fx2(:);      
  res = qdr_intgrl2(lcr%xqdr_cr2DivR2, fx, fx2);
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl_lcr2DivR(lcr, fx, fx2, fx3) result(res);   
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP), intent(IN) :: fx(:), fx2(:);      
  real(DP), intent(IN), optional :: fx3(:);      
  if (present(fx3)) then
    res = qdr_intgrl3(lcr%xqdr_cr2DivR, fx, fx2, fx3);
  else
    res = qdr_intgrl2(lcr%xqdr_cr2DivR, fx, fx2);
  endif
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl_lcr3(lcr, fx, fx2, fx3) result(res);   
  type(GRID_LCR_T), intent(IN) :: lcr;
  real(DP), intent(IN) :: fx(:), fx2(:), fx3(:);      
  res = qdr_intgrl3(lcr%xqdr_cr2, fx, fx2, fx3);
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl(xq, pf) result(res);   
  real(DP), intent(IN) :: xq(:), pf(:);      
!!  call dbg_str("-->qdr_intgrl()", QUADR_TRACE);
  res = vec_dot_prod(xq, pf);
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl2(xq, pf, pf2) result(res);   
  real(DP), intent(IN) :: xq(:), pf(:), pf2(:);      
  !call dbg_str("-->qdr_intgrl2()", QUADR_TRACE);
  res = vec_dot_prod3(xq, pf, pf2);
  return;
end function
!___________________________________________________________
real(DP) function qdr_intgrl3(xq, pf, pf2, pf3) result(res);   
  real(DP), intent(IN) :: xq(:), pf(:), pf2(:), pf3(:);      
  !call dbg_str("-->qdr_intgrl3()", QUADR_TRACE);
  res = vec_dot_prod4(xq, pf, pf2, pf3);
  return;
end function
!___________________________________________________________
! 5-point integration quadratures
! Not tested for dimension(J:J2), where J != 1
SUBROUTINE qdr_step_p5(res, px);      
  real(DP), intent(OUT) :: res(:);      
  real(DP), intent(IN) :: px(:);      
  real(DP) :: stx, h2, wa(4); ! work array; ! step x
  integer nx, mini, maxi, i;
  real(DP) :: INT_P5(4);
  DATA INT_P5 /D14, D32, D12,  D32 /;

  call dbg_str("-->qdr_step_p5()", QUADR_TRACE);
  
  call vec_init_res(res, px);
  call vec_range(mini, maxi, px, nx);
  call dbg_int("mini=", mini, QUADR_TRACE);
  call dbg_int("maxi=", maxi, QUADR_TRACE);
  call dbg_assert_true("xn >= 5", nx >= 5);
  call dbg_assert_true("must be n=k*4+1; MOD(nx-1,4)==0", &
    MOD(nx-1, 4) == 0); ! 
  call grid_step_assert_4(px);
  stx = px(mini+1) - px(mini); ! step
  h2 = D2 * stx / D45;
  wa = h2 * INT_P5; ! load weights for the given x-grid
 
  call vec_dbg("px=", px);
  ! not tested but trying to have independent mini
  do i = mini, maxi
    res(i) = wa(1 + MOD(i - mini, 4)); 
  end do;
  res(mini) = res(mini) * D05;
  res(maxi) = res(maxi) * D05;  
  return;
end SUBROUTINE
!___________________________________________________________
end module 